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1.  Program  Overview 

This  research  program  has  focused  on  the  development  of  novel  numerical  techniques  for  the 
solution  of  electromagnetic  scattering  problems  that  are  accurate,  fast,  and  robust.  The  specific 
advancements  of  this  program  have  been: 

i)  The  derivation  of  a  high-order  method  of  moment  formulation  with  a  quadrature  point 
collocation  scheme  that  is  equivalent  to  a  locally-corrected  Nystrom  (LCN)  method. 

ii)  The  development  of  mixed-order  basis  functions  for  general  LCN  formulations. 

iii)  The  demonstration  of  a  high-order  LCN  formulation  for  the  electromagnetic  scattering  by 
inhomogeneous  material  regions. 

/v)  The  development  of  a  hybrid  volume/surface  integral  equation  formulation  for  the 
electromagnetic  scattering  by  composite  material  objects. 

v)  The  development  of  computer  software  {Mat-Scat)  for  the  simulation  of  the 
electromagnetic  scattering  by  composite  material  objects  based  on  a  high-order  LCN 
formulation. 

vi)  The  a  pplication  o  f  f  ast  iterative  s  olvers  w  ithin  Mat-Scat,  a  nd  t  he  de  velopment  o  f  t  he 
Quadrature-Sampled  Pre-Corrected  FFT  (QS-PCFFT)  formulation. 

vii)  The  validation  of  the  high-order  solution  techniques  for  a  number  of  scattering  problems. 

2.  Accomplishments 

Through  the  course  of  this  program,  the  following  milestones  have  been  accomplished. 

2.A  High-order  method  of  moment  solution  with  Point-Based 
Discretization 

We  have  developed  a  new  method  of  moment  methodology  based  on  high-order  basis 
functions  and  a  quadrature  point  collocation  scheme  [1,2].  By  applying  an  efficient  far-field 
transformation,  the  degrees  of  freedom  can  be  mapped  from  the  basis  function  to  the  quadrature 
points  leading  to  a  form  that  is  identical  to  the  Locally  Corrected  Nystrom  (LCN)  method  [3]. 
The  high-order  method  leads  to  exponential  convergence  for  smooth  scatterers  [2].  Exponential 
convergence  has  also  been  demonstrated  for  geometries  with  edge  singularities  through  the  use 
of  singular  basis  and  specialized  quadrature  rules  [2,  4].  For  arbitrary  geometries,  we  have  also 
demonstrated  that  the  use  of  divergence-conforming  mixed-order  basis  functions  can 
significantly  improve  the  accuracy  of  the  method  and  greatly  improve  the  condition  number  of 
the  impedance  matrix  [5], 
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2.B  High-order  solution  for  the  scattering  by  inhomogeneous  material 
volumes. 

A  method  of  moment  formulation  with  quadrature  point  sampling  for  the  scattering  by 
inhomogeneous  material  volumes  has  been  developed  and  validated.  The  method  has  been 
shown  to  yield  exponential  convergence  for  canonical  and  non-canonical  scatterers  [6-9].  This  is 
the  first  kno  wn  int  egral-equation  b  ased  electromagnetic  s  cattering  s  olution  w  ith  a  hig  h-order 
representation  of  a  material  inhomogeneity. 


2.C  Hybrid  Surface/Volume  integral  equation  formulation 

A  novel  integral  equation  method  has  been  derived  for  a  general  hybrid  surface/volume 
formulation.  The  integral  equation  formulation  allows  for  the  combination  of  arbitrary  material 
surfaces,  including  conducting  as  well  as  penetrable  materials,  and  inhomogeneous  volumes.  A 
computer  program  (identified  as:  Mat-Scat)  implementing  this  technique  using  the  proposed 
high-order  method  of  moment  formulation  has  been  developed.  The  code  has  been  fully 
validated  for  curvilinear  quadrilateral  and  hexahedral  cells  with  orthogonal  mixed-order 
divergence  conforming  basis  functions  and  quadrature  point  collocation. 

2.D  Fast  Solution  Methods 

2.D.1  Quadrature  Sampled  Pre-Corrected  FFT  (QS-PCFFT) 

A  novel  solution  algorithm  referred  to  as  the  “Quadrature  Sampled  Pre-Corrected  FFT”  (QS- 
PCFFT)  algorithm  has  been  developed  [4,  9],  This  method  is  in  the  same  class  of  methods  as  the 
PC-FFT  [10]  and  the  Adaptive  Integral  Method  (AIM)  [11].  The  principal  difference  is  that 
moments  are  never  explicitly  computed.  Rather,  the  projection  from  the  general  discretization  to 
the  uniform  grid  is  performed  efficiently  and  to  high  order  using  Gaussian  quadrature.  The 
method  naturally  preserves  the  high-order  convergence  of  the  method  of  moment  solution,  even 
for  singular  basis  functions  [4],  The  method  scales  as  0(N\ogN)  for  the  scattering  by  material 
volumes  or  planar  structures  and  0(N3/2\ogN)  for  surface  scattering.  The  QS-PCFFT  algorithm 
has  been  implemented  within  Mat-Scat. 

2.D.2  Multilevel  Fast  Multipole  Method 

A  fast  iterative  solver  based  on  the  multi-level  fast  multipole  solution  [12-15]  has  also  been 
implemented  within  Mat-Scat.  The  multi-level  scheme  is  based  on  the  algorithm  proposed  by 
Gyuer  and  Stalzer  [16],  and  the  fast  spherical  filtering  scheme  of  Jakob-Chien  and  Alpert  [17] 
was  implemented  to  accelerate  the  disaggregation  step  of  the  MLFMM  algorithm. 

2.E  High-Order  Mesh  Generation 

The  University  of  Kentucky  scattering  c  ode  Mat-Scat  is  capable  of  handling  surface  and 
volume  elements  of  arbitrary  order.  The  majority  of  commercially  available  mesh  generation 
schemes  are  limited  typically  limited  to  quadratic  or  cubic  elements.  For  the  modeling  of  more 
complex  geometries,  an  automated  tool  must  be  available  for  generating  high-order  meshes. 
Thus,  we  have  developed  a  high-order  mesh  generation  tool.  The  initial  surface  description  and 
low-order  mesh  is  first  entered  into  a  commercial  tool  such  as  SDRC  Ideas™.  The  UK  mesh 
tool  can  convert  the  coarse  low  order  mesh  into  a  high-order  mesh  with  elements  of  arbitrary 
order.  Other  post-processing  tools,  such  as  extrusion  to  an  arbitrary  surface  (e.g.,  for  the 
modeling  o  f  t  hin  m  aterial  c  oatings)  a  nd  m  esh  q  uality  c  hecking  ha  ve  also  b  een  inc  luded.  A 
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graphical  user  interface  (GUI)  has  been  developed  for  the  tool  using  JAVA.  This  makes  the  toll 
convenient  to  use  and  platform  platform  independent. 
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4.  Summary  of  Results 


4.A  High-order  method  of  moment  solution  with  Point-Based  Discretization 
4.A.1  Formulation 


Consider  the  integral  equation  used  to  solve  for  a  surface  current  density  ./  (r') : 


fc{r)=  \K(r,?)J(?)dtf 

s 


(1) 


where  S  is  a  smooth  surface,  <f)mc  (r)  is  the  known  forcing  function  evaluated  at  position  r  on  S, 
K(r,r  )  is  the  kernel.  S  is  discretized  into  Np  curvilinear  patches  that  represent  the  surface 
contour  to  high-order.  The  surface  current  J(r')  is  approximated  over  each  patch  by  a  basis 


function  expansion: 

0  (0  (2> 

A=l 

where,  bk  are  constant  coefficients  and  Fk  (F)  represent  a  set  of  smooth  basis  functions 
distributed  over  the  plh  patch  that  is  complete  to  order  Nk .  A  set  of  smooth  test  functions  T,  (r) 
with  support  over  each  patch  and  complete  to  order  Nk  is  introduced.  The  inner  product  of  (1) 
with  each  of  the  test  functions  is  then  performed,  leading  to: 


lT,(r)r(r)ds  =  \T,(f) 


Nr  Nk 


r=\  k= i 


ds . 


(3) 


On  the  m-th  patch,  the  outer  integral  of  (3)  is  approximated  using  a  Nk  -point  quadrature  rule  that 
integrates  exactly  functions  to  at  least  order  Nk : 


Nk  Nk 

c/=l  I 


NP  Nk 


p-\  k=\ 


(4) 


After  applying  a  simple  linear  transformation  (see  equations  (5)-(7)  in  [1]),  (4)  is  reduced  to: 

r(f„ )  -  ifx  (r')ds’  (5) 

/>=!  <•=!  S 

where  rq  is  a  quadrature  abscissa  point  on  the  mlh  patch.  The  discrete  form  in  (5)  appears  to 

be  a  method  of  moment  formulation  with  a  quadrature  point  collocation  scheme.  However, 
it  is  equivalent  to  a  method  of  moment  procedure  employing  test  functions  up  to  order  Nk 


with  the  outer  integral  evaluated  by  a  fixed-point  quadrature  rule. 

In  the  regions  where  the  patches  are  sufficiently  far  from  the  observation  point,  the  kernel  is 

smooth  and  the  convolutional  integral  in  (5)  can  also  be  approximated  to  high-order  by  a  fixed- 
point  quadrature  rule.  Again,  assuming  an  Nk  -point  rule  that  integrates  smooth  functions 
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exactly  to  at  least  order  Nk ,  the  contribution  to  the  right-hand  side  of  (5)  from  patches  that  are 
sufficiently  far  from  r  can  be  expressed  as: 


L  I*. 

pefar  A'-l 


Nk  (  \  (  \ 
q=\ 


(6) 


where  r  are  the  abscissa  points  and  con  are  the  weights  of  the  quadrature  rule.  In  operator 

Vp  vp 


form,  (6)  can  be  expressed  as: 

_  <7> 

where  the  matrix  [Lp]  is  matrix  local  to  patch  p  with  entries  [Lp]kj  =Fk{rj  ),  the  vector  kqm  is 
the  kernel  K^rqi,rq/i  j  sampled  at  the  quadrature  points ,  and  b  is  the  vector  of  coefficients 


weighting  the  basis  functions.  It  is  then  recognized  that 

<8> 

where  J  are  the  currents  evaluated  at  the  discrete  quadrature  points.  Thus,  (8)  implies  a 
transformation  of  unknowns  from  the  current  basis  to  the  quadrature  points.  This  transformation 
is  then  carried  into  the  near  field  region.  In  this  region,  the  convolutional  integral  in  (5)  must  be 
performed  via  adaptive  quadrature  to  desired  precision.  In  operator  form,  this  is  expressed  as: 

<9> 

where  the  Z-th  element  of  vector  kq  is  =  J  K.(rqn  ,r'\Fk  (r')<Zs',  which  is  evaluated  to  desired 

JSp  m 

precision  using  adaptive  quadrature.  Next,  applying  the  transformation  of  variables  in  (8), 
b  -  [Z,'1  ]7  J  and  (9)  is  rewritten  as: 

<,Kr‘7={[^]X}  7  (10) 

In  summary,  the  contribution  to  the  <7„,-th  row  of  the  impedance  matrix  derived  via  the  method 
of  moment  procedure  for  far  interaction  can  be  expressed  from  (7)  and  (8)  as: 

z;:={k,)t  on 

and  for  near  interaction  from  (10)  as: 


It  is  observed  that  both  (11)  and  (12)  are  identical  to  that  derived  via  the  Locally  Corrected 
Nystrom  (LCN)  method  [2,  3].  Consequently,  the  proposed  high-order  method  of  moment 
(HO-MoM)  procedure  is  equivalent  to  a  LCN  formulation.  However,  it  is  derived  using  a 
more  classical  paradigm  understood  by  most  practitioners. 
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4.A.2  Mixed-Order  Basis  Functions 

When  applying  the  high-order  method  of  moment  (or  LCN)  procedure,  it  is  assumed  that  the 
surfaces  or  volumes  of  the  model  are  discretized  using  general  curvilinear  cells  that  model  the 
surfaces  to  sufficient  accuracy.  The  current  density  is  then  expanded  over  each  cell  with  a  set  of 
local  vector  basis  functions.  As  an  illustration,  we  can  assume  general  curvilinear  quadrilateral 
surface  cells.  A  basis  function  expansion  for  the  current  can  be  represented  in  the  form: 

j,  (y,=o-X;*, =» X),  ov 

.7,  (u' J)  =  b!  l’i2r„  («'R  (/,  -0.J/};*,  -0.J/J),  (14) 

where,  following  the  notation  of  Stratton  [18],  ( u\u 2)  are  the  local  curvilinear  coordinates  of 
the  quadrilateral  cell,  a,  are  the  local  unitary  vectors,  J~g  is  the  Jacobian,  Pt  («)  are  /'-order 
Legendre  polynomials,  and  b['k  are  unknown  constant  coefficients.  The  order  of  the  basis 
functions  are  specified  by  M[ .  Classically,  the  LCN  method  has  been  proposed  using  a  basis 
order  that  is  complete  to  polynomial  order  p  [2,  3].  In  this  situation,  M\  =  M\  —  M\=  M\  =  p . 
This  is  referred  to  here  as  a  polynomial-complete  basis  function  set.  Alternatively,  a  mixed- 
order  basis  function  can  be  posed  that  is  divergence  conforming.  That  is,  the  divergence  of  the 
current  density,  which  is  directly  proportional  to  the  charge  density,  is  complete  to  a  polynomial 
order.  This  is  accomplished  if  the  current  density  is  a  mixed-order.  A  /?-th  order  mixed-order 
basis  function  set  is  then  proposed  with  M\  =M\  =p  +  l  and  M]  -M\=  p . 

The  electric  surface  charge  density  within  the  boundaries  of  a  surface  cell  is: 


Pe  = 


j=--± 

JO>  1?\ 


_1_ 

jo> 


du 2 


(15) 


where  the  range  of  j  and  k  is  governed  by  the  type  of  basis  function  used.  If  a  polynomial 
complete  basis  is  used,  the  charge  density  is  represented  by  an  incomplete  polynomial  expansion. 
In  fact,  the  highest-order  basis  function  pairs  (j/j  =  ku  =  p,  /  =  1,2)  will  produce  surface  charge 
densities  that  are  of  mixed  order.  On  the  other  hand,  the  mixed-order  basis  leads  to  a  surface 
charge  density  that  is  polynomial  complete  to  order  p. 

It  was  shown  in  [19]  by  Qah§kan  and  Peterson  as  well  as  Gedney  [5]  that  the  polynomial- 
complete  basis  can  lead  to  spurious  solutions  for  geometries  with  edge  singularities  and  that  this 
can  be  remedied  through  the  use  of  the  mixed-order  basis.  This  is  somewhat  intuitive  for  the 
following  reason.  If  one  observes  the  charge  density  at  the  singular  edge,  one  expects  the  charge 
density  to  exhibit  the  same  leading-order  singularity  as  the  current  flowing  transverse  to  the  edge 
and  as  the  normal  derivative  of  the  current  flowing  normal  to  the  edge.  This  can  be  well 
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approximated  by  the  mixed-order  basis.  For  the  polynomial  complete  basis  this  is  approximated 
by  the  lower-order  terms,  but  there  are  additional  spurious  charges  from  the  highest-order  basis 
that  will  exist.  For  smooth  scatterers,  this  does  not  occur  and  consequently  either  basis  function 
can  be  used  successfully.  However,  since  the  charge  density  is  over-specified  by  the 

polynomial-complete  basis,  a  poorly  conditioned  system  results  for  the  EFIE  operator. 

When  employing  the  mixed-order  basis  functions  presented,  an  appropriate  testing  scheme 

must  be  implemented  within  the  context  of  the  LCN  method.  There  are  two  choices.  The  first  is 

that  one  can  employ  a  single  quadrature  rule  for  a  cell  and  then  use  a  least  square  solution  for  the 

over  or  under-determined  system  to  calculate  the  local  corrections.  The  other  alternative  is  to 

employ  a  pair  of  mixed-order  quadrature  rules  -  one  for  each  test  vector  -  and  solve  a  square 

matrix  for  the  local  corrections.  We  have  implemented  the  latter  option. 

Consider  the  mixed-order  basis  in  which  the  charge  density  is  polynomial  complete  to  order  p. 
Two  test  vectors  t]  and  t2  are  introduced.  Each  of  the  two  test  vectors  are  associated  with  a 

quadrature  rule  that  is  consistent  with  the  order  of  the  basis  functions  J,  and  J2 .  That  is,  given 

Jt  to  be  of  order  (p  +  l)xp ,  a  (p  +  l)xp  point  quadrature  rule  is  defined  for  the  quadrilateral 

patch  (this  is  done  here  as  a  product  of  one-dimensional  Gauss-quadrature  rules).  The  inner  dot 
product  of  the  integral  operator  with  q  is  then  performed  at  each  discrete  quadrature  abscissa 

point.  Similarly,  the  inner  dot  product  of  the  integral  operator  with  t2  is  performed  at  the 
abscissa  points  of  a  px(p  +  \)  point  quadrature  rule.  In  this  manner,  a  square  linear  system  of 
equations  results  for  the  local  corrections  following  the  procedure  outlined  in  [2,  3]. 


4.B  High-Order  Solution  of  the  Volume  Field  Integral  Equation 

An  inhomogeneous  and  isotropic  dielectric  object  with  relative  permittivity  sr  (r)  contained 


within  a  volume  V  and  situated  in  a  homogeneous  free  space  is  illuminated  by  an  incident 
electromagnetic  field  E'"c  with  an  ej(ot  time  dependence.  The  scattered  field  resulting  from  the 
impinging  wave  can  be  computed  through  the  use  of  equivalent  volume  currents  that  are 
distributed  within  V  and  that  satisfy  the  electric  field  integral  equation  (EFIE)  [20-22]: 


E‘  (a-  )  =  - ~  +  jk0?jo  JJJG(F,  r')  ■  J  (r) dv  ( 1 6) 

jcosn{sr{r)-\)  T 

where  E'  is  the  incident  electric  field  in  the  absence  of  the  dielectric  cylinder,  &()is  the  free 


space  wave  number,  rjo  is  the  free  space  wave  impedance,  G  is  the  dyadic  Green's  function: 


G(r,F)  = 


1  ' 
1  +-TVV 
kl 


g0(r,F) 


(17) 
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gn(r,r')  =  e'jk^  I An\r  —  r'| is  the  free-space  Green's  function,  and  J(r)  is  the  equivalent 

polarization  current  density  distributed  over  V  which  is  defined  as: 

J{r)  =  jcoe0E10'  (r)(e,(r)-l).  (18) 


where  E""  is  the  total  electric  field  in  V. 

The  solution  to  the  EFIE  in  (16)  will  be  obtained  via  the  high-order  method  of  moment 
scheme  p  roposed  in  [2].  T  o  t  his  e  nd,  V is  dis  cretized  w  ith  v  olumetric  curvilinear  cells  t  hat 
represent  the  surface  defining  V  to  high  order.  Within  each  cell,  the  current  is  expanded  using  a 

set  of  smooth  vector  basis  functions  that  are  complete  to  order  p.  A  p-point  Gauss-quadrature 
rule  of  order  2p-l  and  with  abscissas  and  weights  (rm,com)  is  introduced  over  each  volume  cell. 

At  each  quadrature  point  three-independent  test  vectors  are  introduced  tm .  A  point-matched 


formulation  is  obtained  by  sampling  (16)  at  the  quadrature  points,  leading  to: 


l  •£'(')  =  ~ 


-  +  jrioKL  ■  (r')]</v' 


(19) 


jcoe„(sr{rm)-\) 

It  is  shown  in  [2]  that  the  formulation  in  (19)  is  equivalent  to  that  employing  pth  order  test 
functions  and  a  fixed  point  quadrature  rule.  It  is  also  shown  in  [2]  to  be  equivalent  to  a  locally 
corrected  Nystrom  (LCN)  formulation  [3]. 

The  volume  V  defined  by  the  dielectric  scatterer  is  assumed  to  be  discretized  by  curvilinear 
cells.  The  cells  are  currently  assumed  to  be  curvilinear  hexahedron.  Each  cell  is  uniquely 
described  by  a  local  curvilinear  coordinate  system.  The  current  within  each  hexahedron  is 
expanded  via  a  set  of  vector  basis  functions  weighted  by  constant  coefficients  b„  as: 

J(f>tb.UA  <2°) 


where  the  jn  are  the  mixed-order  basis  defined  in  (13)  modified  for  a  three-dimensional 
hexahedral  cell.  The  test  vectors  are  also  chosen  as  the  reciprocal  unitary  vectors  so  that 

tm‘]„  =0  if  m*n. 


A  quadrature  r  ule  is  also  defined  f  or  e  ach  cell.  For  hexahedron,  the  rule  is  conveniently 
defined  via  a  product  rule  of  one-dimensional  Gauss-Legendre  quadratures.  Note  that  a  mixed- 
order  quadrature  rule  is  used  for  each  basis  function,  leading  to  a  square  linear  system  of 
equations. 

To  validate  this  scheme,  consider  the  scattering  by  a  spherical  shell  with  a  continuous  radial 
inhomogeneity.  The  relative  permittivity  of  the  sphere  is  defined  as: 


2.0  + 


(b-af 


7.0,  a<r<b 


(21) 


elsewhere 


where  a  =  0.3  m,  and  b  =  0.5  m.  The  relative  permittivity  thus  varies  from  2.0  at  the  inner  radius 
to  9.0  at  the  outer  radius  with  a  quadratic  profile.  The  shell  was  modeled  with  only  a  single  layer 
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of  curvilinear  hexahedral  cells  (a  total  of  24  cells)  that  exactly  represented  the  inner  and  outer 
surfaces.  The  continuous  dielectric  profile  was  then  modeled  exactly  within  the  curvilinear  cell 
through  the  diagonal  term  in  (19).  The  volume  current  density  and  RCS  were  computed  using 
the  proposed  high-order  scheme.  The  bi-static  RCS  of  the  shell  resulting  from  axially  incident 
vertically  and  horizontally  polarized  waves  is  illustrated  in  Fig.  1  for  two  different  basis  orders 
and  is  compared  against  the  exact  solution.  For  a  basis  order  of3x3x3,  there  is  some  apparent 
error  in  the  RCS  (about  1  dB  on  average).  For  a  basis  order  of4x4x4,  the  predicted  RCS  is 
indistinguishable  from  the  exact  solution  at  the  given  resolution  of  the  graph. 

The  mean  relative  error  in  the  RCS  as  a  function  of  unknowns  for  a  single  layer  of  24 
hexahedral  cells  and  various  orders  is  illustrated  in  Fig.  2.  The  error  is  computed  relative  to  the 
known  analytical  solution.  Superimposed  in  this  graph  is  the  mean  absolute  error  in  the  total 
electric  field,  again  referenced  to  the  analytical  solution.  Note  that  the  maximum  electric  field 
magnitude  in  the  shell  was  about  3  V/m.  Exponential  convergence  is  realized  by  the  high-order 
method. 

The  next  case  studied  is  a  dielectric  ogive  shell.  The  geometry  of  this  object  is  illustrated  in 

Fig.  3.  The  shell  has  a  uniform  thickness  of  0.05  m  and  is  composed  of  a  homogeneous 
dielectric  with£r  =  2.56 .  The  ogive  was  discretized  using  a  single  layer  of  curvilinear 

hexahedron  (as  illustrated  in  Fig.  3).  The  mesh  was  generated  using  a  commercial  mesh 
generator  and  cubic  elements  were  used.  The  bistatic  RCS  of  the  ogive  in  the  6  =  90°  plane  is 
illustrated  in  Fig.  4  for  two  different  discretizations  and  are  compared  to  a  reference  solution  that 
was  finely  discretized.  For  this  simulation,  the  frequency  was  299.796  MHz  and  the  incident 
field  was  incident  on  the  tip  of  the  ogive  (d‘"c  -9Q,f'c  =  0).  Both  vertical  and  horizontal 

polarizations  were  simulated.  The  results  are  indistinguishable  for  the  resolution  presented.  The 
mean  relative  error  in  the  bistatic  RCS  versus  the  total  number  of  unknowns  is  presented  in  Fig. 
5.  For  this  data,  the  error  was  computed  relative  to  a  finely  discretized  reference  solution.  For 
each  curve,  fixed  discretizations  were  used  -  24  cubic  hexahedron  and  56  cube  hexahedron  -  and 
increasing  order.  For  each  data  point,  the  basis  order  is  provided  on  the  graph  (e.g.,  4x4x2 ). 
The  first  two  orders  are  the  transverse  coordinates  and  the  third  is  the  radial.  Overall,  the 
solution  exhibits  exponential  convergence  -  though  there  is  a  stutter  in  the  convergence  for  the 
56-cell  discretization. 
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Fig.  1  Bistatic  RCS  of  the  spherical  shell  with  radial 
inhomogeneous  profile,  for  9mc  =  </>mc  =  0° 
for  V-V  and  H-H  polarizations  in  the  <j>  =  0° 
plane. 


Fig.  2  Mean  Error  in  the  RCS  (relative  error)  and 
the  total  electric  field  (absolute  error)  versus 
the  total  number  of  unknowns  for  a  single 
layer  of  24  hexahedral  shells  and  increasing 
order  for  the  inhomogeneous  spherical  shell. 


Fig.  3  Top  view  of  the  dielectric  ogive  shell  (thickness  =  0.05  m)  and  illustration  of  the  56-hexahedron  mesh. 
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Fig.  4  Bistatic  RCS  of  the  dielectric  ogive  shell  in 
the  6  =  90°  observation  plane,  for 
6int  =  90°,  f,c  =  0°  and  a  frequency  of 
299.796  MHz  for  V-V  and  H-H  polarizations. 
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Fig.  5  Mean  relative  error  in  the  RCS  of  the 
dielectric  ogive  shell  versus  the  total  number 
of  unknowns  for  a  single  layer  of  hexahedral 
shells  and  increasing  order. 
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4.C  Hybrid  Volume/Surface  Field  Integral  Equation  Solution 
4.C.1  Formulation 


Consider  the  electromagnetic  interaction  with  an  inhomogeneous  material  scatterer  made  of  a 
composite  of  penetrable  and  conducting  materials,  as  illustrated  in  Fig.  6.  A  radiating  source 
with  an  e'0”  time  dependence  can  be  placed  either  exterior  or  interior  to  the  inhomogeneous 
object.  Applying  the  equivalence  principal,  equivalent  surface  and  volume  sources  are 
introduced  that  radiate  the  scattered  field.  From  these  sources,  the  integral  equations  will  be 
posed  that  will  allow  for  the  unique  solution  of  the  equivalent  sources. 

Following  Fig.  6,  a  surface  separating  volumes  V.  and  V.  is  denoted  as  Su.  Let  SjL  denote 

the  surface  just  inside  Vh  and  Sf  the  surface  just  inside  Vh  Equivalent  current  densities  are  then 


placed  on  all  surfaces  separating  each  material  volume.  The  surface  current  densities  are  located 
on  either  S* .  or  SF- .  Surface  equivalent  current  densities  on  Sf-  are  posed  as: 


ij 


J:=hxH\  ,  AT  =-fi,xE 


Is*. 


I.s*. 


(22) 


where  »,  is  the  unit  normal  directed  into  Vj .  Equivalent  sources  on  surface  Su  are  defined  as: 


J:  -h.xH 


m : 


■  —Hj  x  E 


(23) 


where  hj  is  the  normal  directed  into  V..  The  upper  and  lower  currents  can  be  constrained  by 
enforcing  the  continuity  of  the  tangential  fields  on  ,  leading  to  the  identity: 

(24) 

where  it  is  recognized  that  n  ■  =  —hr  On  the  surface  of  a  conductor,  only  the  electric  current 


density  is  supported.  Consequently, 


J.  „=h,  xH 


s* 


(25) 


Each  material  volume  is  assigned  a  background  material  profile  (f )h,juih)  .  This  defines  the 
volume  equivalent  currents  within  the  material  volume.  Specifically,  if  (st ,  /u, )  *  ( s,h ,  juih ) 


J V*  —  j®£0£it 


° ir  1 

\Sib  J 


E, 


(26) 


My, 


^-1 

Mib 


H 


where  is  the  relative  permittivity  and  permeability  of  the  physical  material  medium. 

Note  that  if  (el,/ul)  =  (elh,^ilh) ,  the  volume  current  densities  are  identically  zero.  It  is  also  noted 
that  on  any  surface  S(.;,  if  {slh,/uih)  =  (elh,/Jjh J,  no  surface  current  densities  will  be  assigned  to 


that  surface. 
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Fig.  6  Equivalent  currents  representing  the  scattered  fields  due  to  a  plane  wave  impinging  on  an  arbitrary 
inhomogeneous  material  scatterer. 


Each  equivalent  current  density  is  residing  in  an  equivalent  homogeneous  material  volume. 
For  example,  Jf,  M* JK/,  Mv  effectively  radiate  in  a  homogeneous  space  with  material  profile 


(elh,/u,h)  ■  The  equivalent  currents  are  said  to  radiate  a  scattered  field.  The  electric  and  magnetic 


scattered  fields  due  to  the  current  densities  in  V.  can  be  derived  as: 

K"  (l,,K )  =  9,L,  (J„ )-  K,  (M„ )  (27) 

nr  (3*K )  ■ = K.  ( ) + iX  (m„  )  <28) 


where, 


l  )=-;*,/ 


/+-fvv- 

k2 


X"(?)GAr,?)da 


K\X«,)  =  frxGXr,7)Xeq{r')dn' 


(29) 

(30) 


where  Q  is  a  surface  for  a  surface  integral  and  a  volume  for  a  volume  integral,  /  is  the  idem 
factor,  or  unit  dyad,  and  G,  (F,r')  is  the  Green  function  for  a  homogeneous  space  with  material 

profile  (%,/g): 


Contract:  MDA972-01 -1-0022 


13 


University  of  Kentucky 


Final  Report 


DARPA/DSO 


e'(r>r)“to|r-?| 

(31) 

and  k,  =  co^£,hE,h  ■ 

Appropriate  integral  equations  must  be  posed  to  solve  for  the  unknown  equivalent  current 
densities.  Consider  a  material  surface  Su  separating  two  material  regions.  From  (22)  and  (23), 

the  fields  must  satisfy  the  boundary  constraints: 

ft,xE'r  ^-M^-fixE™' 

•Y  i  k'j 

(32) 

-fix  E'f  =  M, .  +h  x  Ef 

I  J  n-  1  J  S' 

1  (./  t.i 

(33) 

h  x  H"“  -J-hx  Hf 

l  /  r>+  /,/  /  /  0  + 

‘Y/  ‘Yy 

(34) 

-h,xH';s_  —  —J,j  +  fi, x  Hjcal 

(35) 

where  E‘"C,H-"C  are  radiated  by  impressed  sources  in  Vi ,  E.c‘",H.cal  are  radiated  by  equivalent 
currents  in  volume  V.  (similarly  in  E),  and  n.  =  -fii  has  been  assumed.  It  is  also  noted,  that  on 

Sf,  the  equivalent  currents  in  Vj  contributing  to  Ef  ,Hf  are  -<7(/  and  -Mi  r 

The  constraints  on  the  interior  and  exterior  boundaries  in  (32)  -  (35)  are  used  to  formulate  the 
integral  equation  for  material  boundaries.  These  equations  can  be  combined  since  only  two 
constraints  are  needed  to  solve  for  the  two  unknowns.  They  must  be  combined  in  a  manner  that 
avoids  spurious  “interior  resonances.”  Two  options  are  available  within  Mat-Scat.  The  first  is 
the  classical  PMCHWT  technique  [23,  24].  This  is  equivalent  to  adding  (32)  with  (33)  and  (34) 


with  (35),  which  leads  to: 

hxEf  —fix  E"’c 

■V,  1  ■ 

h  x  H,nc  -fix  H'f 

'  S\  J  . 


=  -hxE;co' 
=-fixH;cal 


+  fixE'cal  , 

J  -T, 

(36) 

+  fix  H  f  . 

1  s-, 

(37) 

It  is  observed  that  the  diagonal  terms  cancel.  The  diagonal  terms  resulting  from  the  residuals  of 
the  principal  value  integrals  of  the  K  operators  also  cancel.  Thus,  the  PMCHWT  is  a  first-kind 
integral  equation.  Furthermore,  due  to  the  VV  term  of  the  L-operator,  the  PMCHWT  has  a 
hypersingular  kernel.  A  consequence  of  this  is  that  the  operator  is  not  expected  to  be  mesh 
stable.  That  is,  as  the  discretization  is  continually  refined  -  increasing  accuracy  -  the  condition 
number  of  the  discretized  operator  will  grow.  However,  because  of  the  VV  term,  the  PMCHWT 
operator  explicitly  represents  the  quasi-static  field  affects,  and  it  can  be  expected  to  model 
objects  with  geometric  singularities  well. 

We  can  also  look  at  another  interesting  property  of  the  PMCHWT  operator.  Consider  the 
situation  when  lm(f ; )  — >  -jao ,  that  is  the  limit  that  the  material  body  becomes  a  PEC.  It  can  be 


shown  that  in  this  limit,  (36)  will  reduce  identically  to  the  EF1E  and  (37)  will  reduce  to  the  MFIE 
in.  Also,  consider  the  limit  that  — >  oo .  In  this  limit,  the  boundary  becomes  a  hard  electric 
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boundary,  and  the  fields  in  Vj  will  tend  to  0,  and  Mi  j  -^0.  It  can  be  shown  that  in  this  limit, 

the  operators  in  (36)  and  (37)  are  bounded.  Consequently,  the  PMCHWT  operator  is  stable 
relative  to  the  refractive  index  of  the  material  object. 

The  second  option  available  within  Mat-Scat  is  a  Muller  formulation  [25],  This  formulation  is 
derived  by  weighting  (32)  by  sjh,  (33)  by  —ejh  and  adding  the  equations  to  obtain: 

**,»£rL  (38) 


A  _  r?  me 


,  A  „  r^inc 

r  +n,x£jhEj 


r  =-M „(£,„+ eJb)-h, 


*  reseat 

l  ~£J»  >7:XEJ 


Dually,  (34)  is  weighted  by  juih ,  and  (35)  by  -jujh  and  adding  them,  leading  to: 


n,  x n,hHmL I  +  w,  x HjhHmc\  =  J,t,  (//,»  +  Mjh )~n,  x 


S* 


-n,  x 


Hs 


(39) 


The  K-operator  requires  a  principal  value  integral,  which  leads  to  a  well  known  residue  plus 
the  principal  value.  After  performing  the  principal  value  integral,  (38)  and  (39)  are  written  as: 

(40) 

'■v, 

,+  -h,  xMihHSC0'l.  (41) 


", x  £,hE, 


,mcV 


where,  the  surface  integrals  over  use  the  modified  K  integral: 

K  t-X"=  jVxG,(r,F)X„(r')da 


A  _  r"  scat 
r  -  x  £jh  Ej 


(42) 


where  r  e  5  ^.  and  k  =  i  or y. 


It  was  observed  by  Muller  that  for  this  combined  form,  the  hypersingular  static  terms  arising 
from  the  L-operators  will  cancel  with  vanishing  separation  of  the  source  and  observation  points 
([25],  pg.  300).  Thus,  the  effective  singularity  of  the  combined  operator  is  only  1/7?. 
Consequently,  the  Muller  formulation  is  a  second-kind  integral  equation  that  consists  of  a 
diagonal  term  plus  a  compact  operator.  As  a  consequence,  the  Muller  formulation  is  expected  to 
be  well  conditioned  and  mesh  stable  -  that  is  the  condition  number  will  tend  to  a  constant  with 
increased  mesh  discretization.  Since  the  kernel  is  not  hypersingular,  one  can  also  expect  that  for 
smooth  material  bodies,  the  LCN  method  will  converge  exponentially,  and  to  exhibit  a  rate  of 

convergence  s  uperior  t  o  t  he  P  MCHWT  f  ormulation.  A  c  aveat  is  t  hat  in  t  he  ins  tance  w  here 
e .  »  et ,  o  r  /Uj  »  /ui ,  t  he  int  egral  e  quations  in  v  olume  C.  a  re  a  mplified  b  y  t  he  c  onstitutive 


relations.  As  a  consequence,  for  high  contrast  materials,  the  error  in  the  Muller  formulation  is 
amplified. 

In  general,  the  Muller  formulation  is  more  accurate  for  moderate  to  low  contrast  materials 
(er  <  20),  and  the  PMCHWT  formulation  is  more  accurate  for  high-contrast  materials 

Next,  a  on  a  conductor  surface,  one  can  constrain  the  total  tangential  electric  field  to  be  zero 

on  the  surface,  leading  to  the  familiar  electric  field  integral  equation  (EFIE)  [22]: 


ft,  x  E1' 


s* 


= -h,  x  e: 


s* 


(43) 
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The  magnetic  field  is  also  constrained  on  a  PEC  surface  using  (25),  leading  to  the  magnetic  field 
integral  equation  (MFIE)  [22]: 


A  r  rn 

n.  x  H. 


,+  =  Ji,P  ~n,x  H-ca 


s* 


(44) 


The  linear  combination  of  the  EFIE  and  MFIE  is  referred  to  as  the  combined  field  integral 


operator  (CFIE),  which  can  be  expressed  as  [22]: 

«<:/w  *  ■  K"  |s  +  ( 1  -  acm  K 


acnr. 


^  -t-Ef 

) 

■Vo  J 

r  scat 


(45) 


where  aCFIF  is  a  real  constant,  generally  defined  between  0  and  1,  and  t  is  a  vector  tangential  to 
the  PEC  surface  Sip.  When  acm  =  0,  this  reduces  to  the  EFIE,  and  when  aCFlE  =  1,  this 

reduces  to  the  MFIE. 

Next,  in  regions  where  the  volume  currents  are  non-zero,  an  additional  constraint  is  levied. 


EH,K=E-r'(r)  +  Er(r) 

Given  the  volume  currents  defined  in  (26),  (46)  and  (47)  can  be  restated  as: 

jv  (?) 


£“'(?)  = 

'  y  7  FeV, 


i(0£,,£,i 


\eib  j 


H 


T(F)  =■ 

i  V  / 


My 


jm,r,h 


r  r,r 


-  E'cc"  (r  ) 


-Rr(r)- 


(46) 

(47) 

(48) 

(49) 


u*  l) 

Note  that  (48)  and  (49)  are  only  constrained  in  regions  where  Jv  *  0  or  Mv  *  0 .  The  scattered 
field  is  due  to  all  equivalent  currents  present  in  Vi . 


Combining  (36),  (37)  (or  (38)  and  (39)),  (45),  (48),  and  (49)  leads  to  a  highly  versatile  integral 
formulation  f  rom  w  hich  t  he  g  eneral  t  reatment  o  f  t  he  e  lectromagnetic  s  cattering  o  f  a  rbitrarily 
inhomogeneous  material  scatterers  can  be  applied.  The  introduction  of  the  background  material 
within  each  material  volume  can  automatically  eliminate  volume  cells  in  large  homogenous 
regions.  It  also  can  lead  to  a  better  conditioned  formulation  when  simulating  high  contrast 
materials.  This  formulation  is  also  directly  applicable  to  anisotropic  media. 

A  computer  code  based  on  this  integral  formulation  has  been  developed.  The  code  is  referred 
to  as  Mat-Scat  for  general  “material  scattering.”  The  code  has  been  fully  validated  for  surface 
and  volume  scattering  using  curvilinear  quadrilateral  and  hexahedral  cell  discretization.  The 
code  can  handle  an  arbitrary  number  of  material  regions.  Also,  either  the  PMCHWT  or  the 
Muller  formulation  can  be  applied  to  a  material  surface  boundary.  The  software  also  utilizes 
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curvilinear  meshing  of  arbitrary  order.  Mixed-order  basis  functions  are  used  to  expand  the 
current  densities.  While  the  basis  order  can  be  arbitrarily  high,  practically  speaking,  the  highest 
order  that  has  been  simulated  has  been  17th  order.  The  software  is  written  in  FORTRAN,  and  is 
platform  independent. 

4.C.2  Validation 

The  University  of  Kentucky  Mat-Scat  software  has  been  validated  for  a  number  of  canonical 
and  non-canonical  geometries.  In  this  section,  a  representative  set  of  validating  results  are 
presented. 

4.C.2.1  Scattering  by  Material  Spheres 

Initially,  consider  the  scattering  by  a  PEC  sphere  that  is  coated  with  a  thin  dielectric  shell. 
The  PEC  sphere  has  a  radius  a  defined  by  k0a=  3.5.  The  dielectric  coating  has  a  relative 
permittivity  of  er  =4- y'0.1  and  an  outer  radius  b  defined  by  kab  =  3.8 .  The  integral  equation 

solution  employed  Muller’s  formulation  for  the  dielectric  boundary,  and  the  CFIE  with  a  =  0.2 
for  the  PEC  boundary.  Two  meshes  were  used  to  simulate  this  problem.  The  first  used  a  total  of 
48  7,h-order  quadrilateral  cells  to  model  both  the  PEC/dielectric  boundary  and  the  dielectric/air 
boundary.  The  second  mesh  used  a  total  of  108  7th  order  quadrilateral  cells.  For  each  mesh,  the 
order  of  basis  was  increased  from  2  to  6.  Table  1  presents  the  mean  error,  total  CPU  time,  and 

the  condition  number  estimate  for  the  LCN  solution  employing  either  the  polynomial  complete 
or  the  mixed-order  basis  functions  for  the  two  meshes.  In  this  table,  P  is  the  basis  order,  Nc  is 

the  total  number  of  quadrilateral  cells,  and  A  is  the  total  number  of  degrees  of  freedom  (also,  the 
rank  of  the  matrix).  For  all  discretizations  the  mixed-order  basis  exhibits  improved  solution 
error  as  compared  to  the  polynomial  complete  basis.  It  is  observed  that  the  condition  number 
resulting  from  the  polynomial  complete  basis  is  much  larger  than  that  recorded  for  the  mixed- 
order  basis.  This  is  predominately  due  to  the  fact  that  the  polynomial  complete  formulation  is 
not  charge  conserving,  and  the  spurious  charges  in  the  solution  introduce  small  singular  values  in 
the  system  matrix. 


Table  1  Scattering  by  a  dielectric  coated  PEC  sphere  with  PEC  sphere  radius  k()a  -  3.5,  and  outer  shell  radius  k()a  - 
3.8  computed  with  polynomial  complete  (P.  C.)  and  mixed-order  (M,  O.)  basis  functions. 


N 

Mean  Error 

CPU-s 

Condil 

tion  # 

p 

Nc 

P.  C. 

M.O. 

P.  C. 

M.O. 

P.  c. 

M.O. 

P.  C. 

M.O. 

2 

48 

1296 

864 

0.49 

0.29 

98 

69 

123 

45 

4 

48 

3600 

2280 

2. 1x1 0'2 

6.2x1 0'3 

362 

318 

8711 

48 

6 

48 

7056 

6048 

3.4x1 0’3 

3.1x104 

1965 

985 

10377 

51 

2 

108 

2916 

1944 

5.8x1  O'2 

0.10 

231 

190 

126 

43 

3 

108 

5184 

3888 

2.5x1 02 

6.7x1  O'3 

614 

406 

835 

47 

4 

108 

8100 

6480 

3.6x1  O'3 

1.0x1 0‘3 

2392 

979 

6217 

49 

Contract :  MDA972-01-1-0022 


17 


University  of  Kentucky 


Final  Report 


DARPA/DSO 


Another  example  of  a  dielectric  coated  metal  sphere  is  considered.  A  10  cm  radius  PEC 
sphere  is  coated  by  a  1  cm  thick  spherical  dielectric  shell  with  relative  permittivity  er  =  9-y0.3 . 

The  sphere  is  illuminated  by  a  plane  wave  at  1  GHz  (effective  radii  of  k0a  =  2.1,  kt)b  =  2.3 ).  This 

case  was  modeled  using  a  combined  CFIE/SIE  formulation  using  either  the  PMCHWT  SIE 
formulation  or  the  Muller  formulation,  and  CFIE/VIE  formulations.  For  the  CFIE/SIE 
formulation,  a  total  of  48  curvilinear  quadrilateral  cells  were  employed,  24  on  the  PEC  surface, 
and  24  on  the  dielectric  surface.  For  the  CFIE/VIE  formulation,  a  total  of  48  curvilinear  cells 
were  also  used,  24  curvilinear  quadrilaterals  on  the  PEC  surface,  and  24  curvilinear  hexahedron 
were  used  to  model  the  dielectric  shell.  The  error  versus  the  number  of  unknowns  was  recorded 
as  the  basis  order  was  increased.  It  is  observed  that  for  this  case,  the  CFIE/Miiller  formulation 
was  superior  in  terms  of  computational  cost.  However,  the  CFIE/VIE  outperformed  the 
CFIE/PMCHWT  formulation.  This  is  generally  true  for  homogeneous  material  problems  with 
moderately  low  contrast  materials. 

To  study  this  somewhat  further,  consider  the  case  of  the  scattering  by  a  dielectric  sphere.  The 

radius  of  the  sphere  is  10  cm.  The  electromagnetic  scattering  was  predicted  at/  =  2  GHz 
(Aa  ~  15cm,  kna  =  4.2),  when  the  sphere  had  two  different  relative  permittivities:  er  =  1.2  and 

sr  =  36-0.3  j .  The  surface  of  the  sphere  was  discretized  into  54  high-order  curvilinear 

quadrilateral  cells  (8th  order  isoparametric  cells).  The  discretization  was  then  refined  by 

increasing  the  order  of  the  basis  functions  (p-refinement).  The  RMS  error  in  the  RCS  predicted 
by  the  LCN  simulation  versus  the  total  of  unknowns  is  illustrated  in  Fig.  8(a)  for  sr  =1.2  when 

the  basis  order  is  increased  from  2nd-order  through  6th-order.  Similarly,  the  RMS  error  in  the 
RCS  for  sr  =  36-0.3  j  is  illustrated  in  Fig.  8(b),  as  the  basis  order  is  increased  from  2nd  to  7th- 

order.  (Note  that  the  total  number  of  unknowns  for  ;?lh-order  mixed-order  basis  functions  is: 
(/>  + 1) x px4xNc,  where  in  this  example  Nc  =  54  cells  [5]).  For  the  low-contrast  sphere  (Fig. 

8(a)),  both  methods  are  converging  exponentially.  However,  the  Muller  method  is  converging  at 
a  faster  rate  sine  the  PMCHWT  operator  has  a  hypersingular  kernel.  Over  the  whole  range,  the 
Muller  method  is  also  at  a  much  lower  level.  For  the  higher  contrast  dielectric  sphere  (Fig.  8(b)), 
the  Muller  method  again  exhibits  a  faster  rate  of  convergence.  However,  the  level  of  error  is 
higher  for  coarser  discretizations.  Consequently  the  PMCHWT  formulation  reaches  2  digits  of 
accuracy  with  fewer  unknowns  than  the  Muller  formulation.  Due  to  the  rapid  convergence,  the 
Muller  formulation  is  more  efficient  if  additional  digits  of  accuracy  are  needed. 

The  condition  number  of  the  matrices  studied  in  the  examples  of  Figs.  8(a)  and  8  (b)  are 
illustrated  in  Fig.  9.  (The  condition  numbers  reported  were  computed  as  the  ratio  of  the 
maximum  to  minimum  singular  values  of  the  impedance  matrix.)  As  expected,  the  Muller 
method  is  stable  relative  to  discretization  due  to  the  properties  of  the  second-kind  integral 
operator.  The  condition  number  resulting  from  the  PMCHWT  algorithm  grows  linearly  with 
discretization  on  the  log-log  plot.  This  is  due  to  the  hypersingular  kernel  from  the  L  operator. 
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This  s  ame  e  xample  w  as  a  Iso  s  tudied  f  or  a  n  h  -refinement.  T  hat  is ,  t  he  o  rder  o  f  t  he  b  asis 

functions  were  held  fixed  at  third  order  and  the  mesh  was  refined..  The  RMS  error  in  the  RCS  is 
illustrated  in  Fig.  10  for  both  the  cases  when  er  =1.2  and  er  =36-0.3  j  .  On  the  log-log  scale, 

the  Muller  formulation  leads  to  third-order  convergence,  which  is  exactly  what  is  expected  of  the 
third-order  basis  functions.  Optimal  convergence  is  expected  for  this  smooth  geometry  since  the 
operator  is  a  second-kind  operator.  On  the  other  hand,  the  PMCHWT  formulation  converges 
approximately  with  0(h'5).  The  reduction  in  error  convergence  is  due  to  the  hypersingularity 

of  the  kernels,  and  the  additional  derivatives  that  are  effectively  applied  to  the  basis  functions. 

Again,  for  the  low  contrast  sphere,  the  Muller  formulation  is  the  most  efficient.  Whereas  for 
the  high-contrast  sphere,  the  PMCHWT  formulation  is  more  efficient  for  coarser  discretizations, 
and  crosses  over  with  the  Muller  fonnulation  at  3  digits. 


co 
o 
o : 


100  1000  6000 


number  of  unknowns 


Fig.  7  RMS  error  in  the  RCS  of  a  dielectric  coated  PEC  sphere  ( er  =9-  j  0.3 )  and  ktia  =  2.1,  k0b  =  2.3 
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Fig.  8  RMS  error  in  the  bi-static  RCS  for  the  scattering  by  a  dielectric  sphere  versus  the  total  number  of  unknowns 
based  on  calculations  made  via  the  LCN  discretization  of  the  Muller,  PMCHWT,  and  PMCHWT-minus 

(PMCHWT-)  formulations,  (a)  Sr  =  1 .2  (b)  er  =  36  -  0.3  j  . 


Fig.  9  Condition  number  of  the  matrix  arising  from 
the  LCN  discretizationc  of  the  dielectric 
sphere  with  the  PMCF1WT  or  the  Muller 
formulation  versus  the  total  number  of 
unknowns  for  increasing  order  and:  (a) 

8r-\2  (b)  er  =36-0.3  j . 


Fig.  10  RMS  error  of  the  RCS  for  the  dielectric 
sphere  for  fixed  order  (p  =  3)  and  increasing 

number  of  cells  when  (a)  £r=1.2  (b) 

€r  =36-0.3  j . 


4.C.2.2  EMCC  Targets 

A  number  of  the  EMCC  Targets  have  been  simulated  using  Mat-Scat  to  validate  the  accuracy 
and  efficiency  of  the  software.  A  number  of  those  are  reviewed  here. 

Initially,  consider  the  EMCC  wedge  cylinder  plate  (Fig.  1 1  (a)).  The  plate  lies  in  the  z  =  0 
plane  and  the  center  axis  of  the  plate  (through  the  tip)  is  aligned  along  the  x-axis.  The 
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curvilinear  mesh  used  to  represent  the  geometry  is  illustrated  in  Fig.  1 1  (a).  Note  that  curvilinear 
cells  were  only  needed  near  the  curved  circular  boundary.  Linear  quadrilaterals  were  used  to 
model  the  remainder  of  the  target.  The  monostatic  RCS  at  an  elevation  of  1 0°  off  the  z  =  0  plane 
(<9  =  80°)  as  calculated  via  the  LCN  solution  of  the  EFIE  using  mixed-order  basis  functions  is 
presented  in  Fig.  11  (b).  These  results  can  be  compared  to  measured  r  esults  and  computed 
results  published  by  Woo,  et  al.,  in  [26]  and  reproduced  in  Fig.  1 1  (c).  Excellent  agreement  is 
seen,  even  for  the  lowest  order.  The  CPU  time  and  problem  dimensions  are  recorded  in  Table  2. 
The  CPU  times  were  recorded  on  a  single  750  MHz  HP-PA  RISC  processor.  The  total  CPU  time 
includes  the  calculation  of  180  solution  vectors  (1  degree  spacing). 

Next,  consider  the  EMCC  metallic  single  ogive  [27],  The  single  ogive  target  has  a  half  angle 
of  22.62  degrees,  a  total  length  of  10  inches,  and  a  maximum  radius  of  1  inch.  The  major  axis  of 
the  ogive  is  aligned  along  the  x-axis.  The  ogive  approximated  by  72  seventh-order  curvilinear 
quadrilaterals  is  illustrated  in  Fig.  1 2  (a).  The  monostatic  RCS  predicted  using  the  LCN  method 
is  illustrated  in  Fig.  12(b)  for  a  frequency /=  1.18  GHz,  and  in  Fig.  13  (a)  for /=  9  GHz.  In  both 
cases,  the  RCS  was  sampled  in  the  0  =  90°  plane  with  1  degree  increments  in  <j> .  These  results 

can  be  compared  with  measured  data  and  independent  computations  published  in  [27],  and 

reproduced  in  Fig.  12  (c)  for/=  1.18  GHz,  and  Fig.  13  (b)  for /=  9  GHz.  Excellent  agreement  is 
observed.  At/=  9  GHz,  the  ogive  is  approximately  7.6  Xa  in  length  and  0.76  X0  in  radius.  This 

was  modeled  with  288  cells  with  p+ 1=4  and  polynomial  complete  basis.  This  lead  to  a  total  of 
only  9,216  unknowns.  The  total  fill  time  was  839  s  on  a  single  750  MHz  HP  PARISC  processor. 
The  factor  and  solve  time  for  180  angles  was  1,278  s. 

The  next  case  studied  was  the  EMCC  metallic  double  ogive,  illustrated  in  Fig.  14  (a).  The  left 
half  of  the  ogive  has  a  half  angle  of  46.4  degrees  at  the  tip  and  a  half  length  of  2.5  inches.  The 
maximum  radius  of  the  half  ogive  is  1  inch.  The  right  half  of  the  ogive  has  a  half  angle  of  22.62 
degrees  at  the  tip,  a  half  length  of  5  inches,  and  a  maximum  radius  of  1  inch  (where  the  two  half 
ogives  meet).  The  ogive  was  aligned  along  the  x-axis.  The  monostatic  RCS  was  again 
computed  via  the  LCN  method.  A  graph  of  the  RCS  versus  <f>  computed  in  the  0  =  90°  plane  is 
illustrated  in  Fig.  14  (b)  at /=  1.57  GHz,  and  in  Fig.  15  (a)  at /=  9  GHz.  Note  that  ^  =  0 
corresponds  to  the  plane  wave  directed  at  the  22.62  degree  tip.  The  results  are  to  be  compared  to 
measured  data  and  independent  calculations  published  in  [27]  and  reproduced  in  Figs.  14  (c)  for 
1.57  GHz  and  15  (b)  at  9  GHz.  Again,  excellent  agreement  is  realized. 
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(b) 


(c) 


Fig.  1 1  EMCC  wedge  cylinder  plate,  (a)  curvilinear  cell  discretization,  (b)  RCS  computed  via  the  LCN  method  with 
mixed-order  basis  for  p+1  =  5, 6,  and  7,  and,  (c)  results  from  [26]. 


#  cells 

p+1 

N 

Fill 

Total  CPUs 

64 

5 

2,560 

197.3  s 

245.5  s 

64 

6 

3,840 

408.0  s 

535.2  s 

64 

7 

5,376 

696.8  s 

996.0  s 

Table  2.  Dimension  and  Computation  Time  for  the  EMCC  Wedge  Cylinder  Plate  (CPU  time  recorded  on  a  single 
750  MHz  HP-PARISC  processor) 


Next,  the  RCS  of  the  metallic  double  ogive  coated  by  a  thin  lossy  dielectric  coating  was 
analyzed.  The  geometry  is  illustrated  in  Fig.  16  (a).  The  coating  was  added  on  by  extending  the 
surface  of  the  double  ogive  by  10  %  along  the  surface  normal.  Thus,  the  maximum  radius  of  the 
coating  is  1.1  inches.  The  left  half  length  of  the  coated  ogive  is  2.75  inches,  and  the  right  half 

length  of  the  coated  ogive  is  5.5  inches.  The  coating  was  assumed  to  be  isotropic  homogenous 
with  a  relative  permittivity  er  =9-y'0.3 .  The  monostatic  RCS  computed  in  the  6  =  90°  plane 
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is  illustrated  in  Fig.  16  (b)  as  computed  via  the  LCN  with  a  combined  field  surface  integral 
equation  formulation,  as  defined  in  Section  2.1.  For  this  discretization,  the  RCS  has  converged 
to  <  1  %  accuracy  with  fourth  order  basis. 

Next,  the  RCS  of  the  NASA  Almond  is  studied  [27],  The  NASA  almond  is  a  challenging 
structure  in  that  the  geometry  is  singular  at  the  two  tips  and  there  is  a  high-rate  of  curvature  in 
the  geometry.  An  example  of  a  curvilinear  mesh  used  to  model  the  almond  is  illustrated  in 
Fig.  17(a).  The  monostatic  RCS  of  the  almond  was  computed  in  the  zero  elevation  plane 
(0  =  90°)  as  a  function  of  the  azimuthal  angle  <f> .  The  broad  side  of  the  almond  lies  in  this 
plane,  and  <j)  =  0°  corresponds  to  incidence  on  the  sharp  tip.  The  RCS  of  the  Almond  predicted 
by  the  LCN  method  at  7  GHz.  The  LCN  simulation  employed  a  140-cell  eighth-order 
curvilinear  mesh  and  mixed-order  basis.  The  RCS  predicted  by  Mat-Scat  is  illustrated  in  Fig.  17 
(b)  for  horizontal  polarization  (HH)  for  fifth  and  sixth-order  discretizations,  requiring  5,600  and 
8,400  unknowns,  respectively.  The  difference  between  these  two  simulations  is  <  .001  dB.  A 
Galerkin  simulation  employing  low-order  basis  functions  was  also  used  to  simulate  the  Almond. 
The  Galerkin  simulation  utilized  a  mesh  consisting  of  21,120  bi-linear  quadrilaterals,  supporting 
31,680  zeroth-order  divergence  conforming  GWP0  basis  functions.  The  low-order  Galerkin 
simulation  still  has  not  converged  for  angles  of  incidence  near  the  tips  of  the  almond,  whereas 
the  LCN  solution  has  converged  with  far  fewer  unknowns.  Figure  17(c)  compares  the  RCS  as 
predicted  by  the  LCN  solution  of  the  CFIE  with  mixed-order  basis  versus  that  computed  with 
polynomial-complete  basis  functions.  Both  simulations  used  the  same  discretization  as  that  in 
Fig.  7  and  sixth-order  basis  functions.  The  polynomial  complete  basis  still  results  in  significant 
error  for  incidence  near  the  tips.  These  results  can  be  compared  against  measured  data  and 
independent  computations  published  by  Woo  et  al.  in  [27],  and  reproduced  in  Fig.  17  (d). 

Finally,  consider  the  case  of  the  thick  trapezoidal  plate  [28,  29]  coated  along  the  edges  with  a 
thin  dielectric  coating  (see  Fig.  18).  This  geometry  has  also  been  referred  to  as  the  EMCC  Dart. 
The  base  of  the  PEC  plate  is  7  feet  long  and  3  feet  high.  The  top  of  the  plate  is  3  feet  long.  The 

plate  is  1  inch  thick.  The  edge  of  the  plate  is  coated  with  uniform  dielectric  coating  with  relative 
permittivity  sr  =  4.5-  j9 .  The  coating  is  2  inches  wide  (except  at  the  comers)  and  is  1  inch 

thick.  The  plate  and  coating  were  modeled  with  a  total  of  158  bi-linear  quadrilaterals.  The  plate 
is  situated  in  the  z  =  0  plane.  The  monostatic  RCS  of  the  plate  was  computed  via  the  LCN 
method  with  mixed-order  basis  in  the  zero  degree  elevation  plane  (0  =  90° ).  The  angle  <f  =  0° 

represents  incidence  on  the  sharp  tip  of  the  plate.  The  RCS  is  illustrated  in  Figs.  19  (a)  and  (b). 
These  results  can  be  compared  to  Fig.  7  of  [29].  The  vertical  polarization  converges  quite  well 
for  order  q  =  5  (N  =  9,200).  Whereas,  the  horizontal  polarization  requires  q  =  6  for  convergence 
(N=  13,800). 
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(a) 


Fig.  12  EMCC  metallic  single  ogive,  f  =  1.18  GHz.  (a)  Patch  discretization  -  72  seventh-order  curvilinear  cells,  (b) 
RCS  computed  via  the  LCN  method  for  p+ 1  =3, 4,  and  5,  (c)  results  from  [27], 


Azimuth 


(a) 


(b) 


Fig.  13  EMCC  metallic  single  ogive,  /=  9  GHz,  (a)  RCS  computed  via  the  LCN  method  for  p+ 1  =4,  and  288 
seventh-order  curvilinear  cells,  (b)  results  from  [27]. 
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(b)  (c) 

Fig.  14  EMCC  metallic  double  ogive, /=  1.57  GHz,  (a)  discretization  with  56  seventh-order  curvilinear  cells,  (b) 
RCS  computed  via  the  LCN  method  for  p+ 1  =4,  and,  (c)  results  from  [27]. 


-10 


A  .  „  0  2i>  10  fiu  W  100  I2<>  NO  N>*.j  180 

Azimuth  A/imuih 


Fig.  15  EMCC  metallic  double  ogive,  f  =  9  GHz,  (a)  RCS  computed  via  the  LCN  method  with  224  cells  and  p+1  = 
4,  (b)  results  from  [27]. 
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(a) 


(b) 


Fig.  16  Coated  EMCC  metallic  double  ogive, /=  1.57  GHz.  Coated  with  thin  dielectric  coating  with  sr  =  9-y'0.3  . 
(a)  Curvilinear  discretization,  (b)  RCS  (monostatic)  computed  in  the  6  =  90°  plane  via  the  LCN  method  for 
p+\  =  4,  5,  and  6. 


(a) 


♦ 


♦ 

(b) 


.1*- 


UC  UO  if*  ;,40 


(C)  (d) 

Fig.  17  NASA  Almond,  /=  7  GHz,  (a)  discretization,  (b)  RCS  computed  via  the  LCN  method  for  p+ 1=3  via  the 
CFIE  and  MFIE  with  polynomial  complete  basis  as  well  as  a  Galerkin  scheme  with  RWG  basis  (HH),  and, 
(c)  LCN  method  with  p  =  5  via  the  EF1E  with  mixed-order  basis  as  compared  to  a  Galerkin  method  with 
RWG  basis  (HH),  (d)  results  from  [27], 
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Fig.  18  Coated  trapezoidal  plate  -  0.0254  m  thick  PEC  plate  with  a  0.0508  m  dielectric  edge  coating  er  =  4.5  -  j9  . 
The  plate  is  located  in  the  z  -  0  plane. 


<t>  4> 

(a)  (b) 

Fig.  19  Monostatic  RCS  of  the  coated  trapezoidal  plate  computed  via  the  LCN  method  in  the  zero  elevation  plane  at 
1  GHz.  (a)  V-V  polarization,  (b)  H-H  polarization. 
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4.D  Fast  Solution  Methods 

4.D.1  The  Quadrature  Sampled  Pre-Corrected  FFT  Algorithm 


We  have  developed  a  novel  fast,  high-order  solution  procedure  referred  to  as  the  Quadrature 
Sampled  Pre-Corrected  FFT  (QS-PCFFT).  The  method  accelerates  far-interaction  terms  of  an 
integral  operator  using  the  discontinuous  FFT  [30],  which  combines  Gaussian-quadrature 
integration  with  the  FFT.  We  have  applied  the  method  to  our  high-order  method  of  moment 
scheme  with  point-based  discretization  (equivalent  to  the  locally  corrected  Nystrom)  for 
electromagnetic  scattering.  This  is  briefly  described  in  the  following. 

Consider  the  electromagnetic  scattering  by  a  PEC  object.  The  electric  field  integral  equation 
(EFIE)  can  be  used  to  solve  for  the  currents  induced  on  the  PEC  surface  S.  The  EFIE  is 
described  as: 

E'(r)=  Jp(r,r')J(r')</r'  (50) 


where  E'  is  the  incident  electric  field,  ko  and  rj0  are  the  free  space  wave  number  and  wave 
impedance,  G  is  the  dyadic  Green's  function: 


G(r,r')  =  j?]X 


I+- yVV 
kl 


go(r>  r') 


(51) 


g0(r,r')  =  e  Jk^r  / 4zr  r-r'  is  the  free-space  Green's  function,  and  J(r)  is  the  equivalent 


current  density.  Using  the  method  of  moment  procedure  detailed  in  [2]  leads  to  a  discrete  linear 
system  of  equations: 

e  =  Zb  (52) 

where  Z  is  a  dense  matrix,  e  represents  the  known  forcing  vector  and  b  is  the  unknown  solution 
vector.. 

It  is  assumed  that  the  linear  system  in  (52)  will  be  solved  using  an  iterative  scheme.  Thus, 
each  iteration  will  require  a  matrix  vector  product.  This  will  be  accelerated  using  the  proposed 

fast  scheme.  Initially,  the  impedance  matrix  is  decomposed  into  two  parts: 

Z  =  Zncar  +ZM .  (53) 

where,  Znear  represent  the  near-field  blocks  that  are  computed  to  desired  precision  using 
adaptive  quadrature.  Zfar  represents  the  far-field  blocks  that  are  computed  efficiently  using 
single-point  evaluations  of  the  kernel.  The  proposed  fast  scheme  will  accelerate  the  product  of 
the  matrix-vector  product  involving  Z  far .  This  is  described  in  the  following. 

Due  to  the  translational  invariance  of  the  Green’s  function,  the  convolution  of  the  current  with 
the  d  yadic  G  reen’s  f  unction  c  an  b  e  computed  e  fficiently  f or  far  int  eractions  v  ia  t  he  dis  crete 
Fourier  transform  (DFT).  Here,  the  DFT  is  performed  using  the  “discontinuous-FFT”  (D-FFT) 
recently  introduced  by  Fan  and  Liu  [30].  This  method  allows  one  to  compute  the  FFT  of  a 
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discontinuous  and/or  singular  function  to  controllable  accuracy.  Since  the  method  is  based  on  a 
quadrature  sampling,  it  is  ideal  for  a  Nystrom  discretization. 

Fan  and  Liu’s  discontinuous  FFT  method  is  extended  here  to  multiple  dimensions.  Consider 
the  FFT  of  a  multi-dimensional  function  /( r)  distributed  on  Q.  Q.  is  discretized  into  L  general 

curvilinear  cells,  as  governed  by  the  method  of  moment  discretization.  The  multi-dimensional 
discrete  Fourier  transform  of  /( r)  can  be  defined  as: 

/(")  =  !  f  /(r)e«“Vn  (54) 

/=i  n7 

where  n  =(«,,«2,...)  defines  the  index  of  the  discrete  Fourier  domain.  The  integrals  are  then 


approximated  via  an  appropriate  numerical  quadrature  rule: 


l=\ 


(55) 


where  is  the  Jacobian  evaluated  at  the  abscissa  point  r'q  on  the  /- th  curvilinear  cell,  and 
fq  =  f(r'q) .  A  uniform  grid  with  spacing  Ax  is  superimposed  over  the  general  discretization  (see 


Fig.  7).  The  indices  of  the  multi-dimensional  grid  is  defined  to  have  physical  coordinates 
rm  =(m,Ax,  m2Ay,...).  The  exponential  function  at  the  indices  of  the  uniform  grid  is  expressed 

as:  eJ2’Tnr"' .  Since  the  exponential  function  in  (55)  is  a  smooth  function  over  all  space,  it  can  be 
interpolated  from  the  uniform  grid  indices  to  the  quadrature  points  as: 

,  M 

e;2'nr,=Sa (56) 


where  <t>m  are  smooth  interpolation  polynomials.  (It  is  noted  that  the  summation  in  (56)  is 
taken  over  the  entire  uniform  grid.  However,  the  interpolation  polynomials  d>)n(r^  are  non¬ 
zero  only  in  a  region  local  to  r‘ ).  Then,  combining  (56)  with  (55),  leads  to 


/=!  q=\  m=0 


ej2™r" 


(57) 


The  order  of  summation  is  then  rearranged  as: 


M  I  ", ,  - - 

m = )f. 


m- 0 
M 


/= 1  ij=\ 


(58) 


=  Y,g,e,M'- 


m= 0 


Consequently,  /(n)  is  evaluated  efficiently  for  all  values  of  n  in  the  discrete  space  via  the  FFT. 


The  discontinuous  FFT  in  (58)  is  used  to  evaluate  a  matrix-vector  product  involving  Z  far .  Let 
j  be  the  vector  of  currents  at  the  quadrature  points.  Then,  the  FFT  of  j  is  expressed  in  operator 
format  as 


j(n)  =  «T{0]}. 


(59) 
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where  the  operator  W  is  a  highly  sparse  matrix  with  elements  governed  by  the  double  summation 
in  (58)  and  is  the  multi-dimensional  FFT.  Next,  the  scattered  field  due  to  j  can  be  computed 
at  the  uniform  grid  points  using: 

=  -HW] 

where  the  dyadic  Green’s  function  G  is  evaluated  at  the  uniform  grid  points  referenced  to  the 
origin.  The  field  must  be  evaluated  at  the  quadrature  points.  Since  the  field  away  from  the 
source  is  smooth,  it  can  be  interpolated  using  the  polynomials  j .  Consequently, 

e*(r 'q)  =  -VTHW\.  (61) 

where  IF  is  a  sparse  matrix.  Furthermore,  it  can  be  shown  that  W  =  ci)V ,  where  a>  is  a  diagonal 
matrix  with  the  <7/-th  entry  . 


Finally,  we  can  express: 


ZMj  =  VtHW] 


(62) 


where  j  is  the  vector  of  currents  at  the  quadrature  points. 

Two  items  should  be  realized  at  this  point:  1)  the  far  field  computation  is  limited  to  regions 
where  the  fixed  point  quadrature  provides  sufficient  accuracy  to  the  convolutional  operator,  and 
2)  the  dyadic  Green’s  function  is  hypersingular  at  vanishing  separation  of  the  source-observation 
point  pair.  Since  the  FFT  is  only  applied  in  the  far-region,  it  is  sufficient  to  zero  out  the  dyadic 
Green’s  function  in  the  vicinity  of  the  origin.  Next,  since  (62)  is  a  global  operator  it  will  lead  to 
non-zero  fields  in  the  “near-region”.  This  error  is  compensated  for  by  “pre-correcting”  Znear 
[10]: 

Z"a,r  =  Znear  -  VT HW  (63) 

The  system  matrix  in  (53)  is  then  expressed  as: 

Z  =  Z"car  +VTHW .  (64) 

This  fast  scheme  is  referred  to  here  as  the  QS-PCFFT  method.  This  technique  can  be 
compared  to  both  the  PC-FFT  [10,  31,  32]  and  AIM  methods  [11,  33-35],  The  most  distinctive 
difference  between  these  techniques  and  the  QS-PCFFT  is  the  projection  to  the  uniform  grid. 
Both  the  AIM  and  the  P  C-FFT  methods  equate  “moments”  computed  via  monomials  on  the 
uniform  grid  and  the  current  basis.  This  computation  is  somewhat  analogous  to  computing  the 
weights  of  a  numerical  quadrature  rule.  However,  the  Gaussian  quadrature  rule  employed  by  the 
D-FFT  inherently  integrates  these  moments  to  higher  order  -  even  for  singular  or  discontinuous 
functions  using  the  appropriate  quadrature  rule.  In  fact,  the  weights  and  abscissa  are  chosen  to 
realize  optimal  convergence. 

It  is  next  recognized  that  the  QS-PCFFT  computes  the  far  interactions  to  the  same  precision  as 
the  fixed  point-quadrature  integration.  Thus,  the  precision  of  the  far-field  calculation  is  known 
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relative  to  a  more  rigorous  adaptive  quadrature  evaluation.  The  break  between  the  near  and  far 
interaction  is  thus  specified  by  the  quadrature  order  and  a  desired  error  tolerance. 

The  complexity  and  memory  requirements  of  the  QS-PCFFT  will  scale  in  a  similar  manner  as 
the  PC-FFT  and  AIM  methods.  Since  Z"ear  is  bound  near  the  diagonal,  filling  and  storing  this 
matrix  will  scale  with  problem  dimension  as  O(N)  (for  a  fixed  order).  Furthermore,  pre¬ 
correcting  the  near  impedance  matrix  will  also  scale  as  0(N )  since  it  is  a  local  operation. 

Iterative  solutions  require  a  matrix  vector  multiply,  which  is  performed  with  the  operator  in 
(64).  The  product  with  Z"car  will  scale  as  O(N) .  The  scaling  of  the  product  performed  via  the 

operator  VT HW  will  be  dependent  on  the  distribution  of  quadrature  points.  If  the  quadrature 
points  densely  populate  the  discrete  space  so  that  the  uniform  FFT  grid  dimension  scales  linearly 
with  the  number  of  unknowns,  then  VT HW  will  have  a  complexity  that  scales  as  0(N  log  N). 

However,  if  the  quadrature  points  sparsely  populate  the  discrete  space  (e.g.,  a  three-dimensional 
metallic  surface  such  as  a  sphere),  then  the  complexity  can  scale  as  0{Ny  log  N) .  Thus,  the  QS- 

PCFFT  method  is  most  economically  applied  to  problems  that  involve  a  more  densely 
distribution  of  unknowns  -  that  is,  where  the  unknowns  occupy  much  of  the  problem  space. 

The  QS-PCFFT  is  an  approximation  of  the  convolution  rather  than  an  asymptotic  expansion  of 
the  kernel.  Therefore,  the  dimension  of  Znear  is  not  a  function  of  electrical  lengths,  but  rather  a 
function  of  discretization.  Consequently,  the  method  can  be  applied  to  both  low  and  high 
frequency  problems  with  equal  effectiveness. 

To  demonstrate  the  efficiency  of  the  proposed  scheme,  consider  the  electromagnetic  scattering 
by  an  array  of  circular  rings.  The  unit  cell  of  this  array  (illustrated  in  Fig.  8)  is  situated  in  a  z  = 
constant  plane  and  has  a  dimension  of  1.5  X  x  1.5  X.  The  curvilinear  cells  used  to  represent  the 
rings  are  illustrated  in  Fig.  9.  The  interpolating  grid  of  a  6th  order  polynomial  cell  is  also 
illustrated  for  one  of  the  cells  in  ring  #2.  Jacobi  polynomials  with  singular  weights  that 
appropriately  model  edge  singularities  were  used  for  basis  functions.  A  Gauss-Jacobi  quadrature 
rule  was  used  for  the  Nystrom  formulation  and  the  QS-PCFFT.  Since  the  quadrature  rule 
integrates  the  singular  basis  to  high  order,  the  QS-PCFFT  is  still  exponentially  convergent 
-  even  for  singular  basis  functions. 

The  bi-static  RCS  of  a  10  x  10  array  of  the  unit  cells  is  illustrated  in  Fig.  10  computed  with  the 
discretization  in  Fig.  9,  with  3rd  order  basis  and  3x3  point  quadrature  rules  on  each  patch.  A 
512  x  512  point  FFT  was  used  for  the  QS-PCFFT  solution,  with  the  uniform  grid  being  three 
times  the  global  dimension  of  the  array  (45  X  x  45  X).  The  error  is  estimated  to  be  ~1  %  for  this 
discretization.  The  RCS  was  also  computed  using  FISC  [36]  employing  RWG  basis  on  linear 
quadrilateral  cells  [13,  37].  Two  different  discretizations  were  used  for  the  FISC  simulations 
resulting  in  17,100  (~  10  cells  per  wavelength)  and  113,400  unknowns  (~25  cells  per 
wavelength).  The  results  of  the  higher  discretization  is  much  closer  to  the  LCN  solution. 
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The  CPU  times  and  memory  of  the  simulations  on  a  single  750  MHz  HP-PARISC  processor 
are  recorded  in  Table  3.  (Note  that  the  LCN/QS-PCFFT  simulation  used  double  precision,  and 
the  FISC  simulation  used  single  precision.)  It  is  observed  that  as  the  discretization  of  the  FISC 
simulation  was  increased,  the  memory  increased  more  dramatically.  The  reason  for  this  is  that 
since  the  finest  level  of  the  FMM  is  fixed.  Consequently,  the  number  of  cells  in  the  finest  groups 
increases.  This  increases  the  dimension  of  the  block  matrix  that  must  be  stored. 

To  study  the  scaling  of  the  complexity  of  the  QS-PCFFT  solution,  the  dimension  of  the  array 
of  circular  rings  in  Fig.  8  was  increased  from  lxl  to  20x20.  The  discretization  used  for  each 
unit  cell  was  fixed  to  3rd  order  basis  and  3x3  point  Gauss-Jacobi  quadrature  rules.  The  size  of 
the  FFT  was  also  scaled  with  the  problem  domain.  The  CPU  times  required  for  the  matrix  fill 
and  a  matrix-vector  multiplication  are  presented  in  Fig.  10.  Also  overlaying  these  curves  are 
curves  illustrating  0(N)  and  0(NlogN)  complexity  of  the  tasks,  respectively. 


Fig.  20 


Ax{ 


Uniform  grid  with  spacing  Ax  indices  x/j(  overlaying  curvilinear  cells  discretized  by  quadrature 
abscissa  points  . 


Fig.  21  Illustration  of  the  curvilinear  cells  used  to  model  one  unit  cell  of  the  ring  array.  The  center  coordinates 
(X > X- )  /  Inner/Outer  Radius  (r,,r0)  of  the  five  rings  are:  #1  (0.52, 0.52)/(0. 12,0.252) , 
#2  (0.5/1, 0.5/1) /(0. 35/1, 0.5/1) ,  #3  (0.22, 1.22)/(0.12, 0.22) ,  #4  (1.22,1 .22)/(0.12, 0.22) ,  #5  (1.22,0.22) 
/(0. 12, 0.22)  One  patch  on  ring  #2  illustrates  the  discretization  of  a  sixth-order  polynomial  cell. 


Contract :  MDA972-01-1-0022 


32 


University  of  Kentucky 


Final  Report 


DARPA/DSO 


Fig.  22  Bistatic  RCSofthe  ring  array  with  10x10  unit  cells  (a  =  b  =  1.5  X)  with  (0',^')  =  (0°,0°)  computed  via  the 
present  method  as  compared  and  using  FISC  with  RWG  basis  and  the  MLFMA  [36],  (a)  V-V  polarization, 
(b)  H-H  polarization. 


Fig.  23  CPU  time  for  the  matrix  fill  and  the  time  per  BICGSTAB  iteration  versus  the  number  of  unknowns  recorded 
on  a  single  750  MHz  HP-PARISC  processor  as  the  ring  array  is  scaled  from  lxl  to  20x20  array.  Also 
overlaid  are  curves  with  slope  A  and  slope  of  N\ogN. 


Unknowns 

gillie 

Solve  Time 

Memory 

Total  CPU 

QSPCFFT 

43,200 

658  s 

415s 

109  MB 

1,074  s 

FISC 

17,100 

124  s 

1,562  s 

32  MB 

1,686  s 

FISC 

113,400 

7,930  s 

7,930  s 

331  MB 

11,853  s 

Table  3.  Scattering  by  a  10  x  10  circular  ring  array  recorded  on  a  single  750  MHz  HP-PARISC  processor 
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4.D.2  The  Multilevel  Fast  Multipole  Method  (MLFMM) 

A  fast  iterative  solver  based  on  the  multi-level  fast  multipole  solution  [12-15]  has  also  been 
implemented  within  Mat-Scat.  The  MLFMM  offers  a  more  general  fast  iterative  solution 
algorithm  that  will  scale  as  0(N\ogN)  for  general  three-dimensional  scattering  problems.  The 

multi-level  scheme  implemented  is  based  on  the  algorithm  proposed  by  Gyuer  and  Stalzer  [16], 
and  the  fast  spherical  filtering  scheme  of  Jakob-Chien  and  Alpert  [17]  was  implemented  to 
accelerate  the  disaggregation  step  of  the  MLFMM  algorithm.  For  multi-region  problems,  a 
separate  multipole  expansion  is  performed  in  each  background  region. 

As  an  illustration,  consider  again  the  EMCC  metallic  ogive  in  Fig.  12(a).  Both  the  MLFMM 
and  the  QS-PCFFT  fast  iterative  solvers  were  used  to  simulate  the  RCS  at  9GHz.  The  surface 
was  modeled  with  288  curvilinear  surface  cells.  The  surface  currents  were  modeled  with  5th 

order  basis  -  leading  to  a  total  of  1 1,520  unknowns.  For  the  MLFMM,  the  smallest  group  size 
was  set  to  0.5  A^.  The  QS-PCFFT  used  a  uniform  FFT  grid  that  had  a  dimension  of 

1 08  x  20  x  20 .  The  CPU  times  and  memory  required  by  the  fast  algorithms  are  compared  with  a 
direct  solve  in  Table  4.  A  small  enough  problem  was  solved  to  compare  against  a  direct  solution 
algorithm.  It  is  also  noted  that  the  results  simulated  via  fast  solvers  is  as  accurate  as  that  solved 
via  the  direct  solver.  That  is,  both  the  MLFMM  and  the  QS-PCFFT  preserves  the  accuracy  of 
the  simulation. 


Filling  time  (S) 

Solving  Time  (S) 

Memory  Usage  (Mb) 

Direct  Solver  (LU) 

1,049 

2,239 

2025 

QS-PCFFT 

404 

450 

46.7 

MLFMA 

315 

540 

71.50 

CPU  time  and  Memory  usage  comparison  of  fast  algorithms  with  direct  solver  at  9GHz. 
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4.E  High-Order  Mesh  Generation 

High-order  methods  such  as  the  LCN  method  are  most  efficient  when  employing  higher-order 
basis  o  n  1  arge  s  mooth  c  urvilinear  p  atches.  T  he  r  eason  f  or  t  his  is  s  imple:  hig  her-order  b  asis 
converge  more  rapidly  than  lower  order  basis.  Thus,  with  the  LCN  method  it  is  desirable  to 
model  geometries  with  large  curvilinear  cells  that  represent  the  surface  to  sufficient  accuracy.  In 
practice,  cells  with  average  edge  lengths  on  the  order  of  one  wavelength  are  quite  typical.  Thus, 
it  becomes  imperative  to  employ  isoparametric  curvilinear  patches  that  accurately  model  a 
surface  of  arbitrary  curvature.  If  one  can  not  support  such  patches,  then  one  is  forced  to  use  a 
refined  discretization  -  thus  losing  the  advantage  of  the  high-order  method.  For  classical  low- 
order  techniques,  this  has  not  been  an  issue  since  the  slow  convergence  requires  one  to  resolve 
the  surface  with  a  minimum  of  10  to  20  edges  per  linear  wavelength  to  get  reasonable  accuracy. 
Often,  such  fine  sampling  is  also  enough  to  represent  a  curved  surface  to  sufficient  accuracy 
using  a  piecewise  linear  approximation. 

Most  commercial  mesh  generation  programs  are  limited  by  the  order  of  curvilinear  cells  that 
can  be  generated.  Most  CAD  packages  provide  linear  triangles  and/or  bi-linear  quadrilaterals. 
Some  will  generate  bi-quadratic  and  very  few  will  generate  up  to  bi-cubic  triangle  or 
quadrilateral  elements.  Unfortunately,  this  is  insufficient  for  a  true  high-order  method,  which 

requires  discretization  to  arbitrary  order.  This  is  illustrated  through  a  simple  example.  Consider 
the  electromagnetic  scattering  by  a  PEC  sphere  of  radius  a  defined  by  koa  =  6  .  The  sphere  was 

discretized  with  24  quadrilateral  curvilinear  cells  as  illustrated  in  Fig.  24(a).  The  basis  function 
order  was  fixed  at  p  =  9,  leading  to  3,456  unknowns.  The  bi-static  RCS  was  then  predicted  by 
Mat-Scat  as  the  order  of  the  curvilinear  cell  was  increased  from  n  =  1  to  n  =  12.  A  graph  of  the 
error  in  the  RCS  relative  to  an  exact  Mie-series  solution  is  illustrated  in  Fig.  24(b).  A  very 
important  observation  is  that  despite  using  high-order  basis  functions  (9lh-order),  the  error  in  the 
RCS  is  limited  by  the  error  in  the  discretization  of  the  geometry.  That  is,  if  we  were  limited  to 
bi-cubic  quadrilateral  cells  (3rd  order),  then  the  simulation  would  only  provide  2  digits  of 
accuracy  despite  using  9th-order  basis  functions.  Whereas,  minimum  error  is  realized  with  a 
mesh  that  is  of  order  n  =  p  + 1  (10th  order),  resulting  in  6  digits  of  accuracy. 

The  error  in  the  area  of  the  sphere  as  approximated  via  a  9x  9 -point  Gauss-Legendre 
quadrature  rule  on  each  patch  is  also  graphed  in  Fig.  24  (b).  It  is  very  interesting  to  observe  that 
the  errors  in  the  RCS  and  the  area  follow  the  same  general  trend.  This  shows  that  there  is  a 
direct  correlation  in  the  accuracy  of  the  surface  model  and  the  accuracy  of  the  EM  scattering  as 
predicted  by  the  LCN  method.  Thus,  to  accurately  and  efficiently  predict  the  scattering  by  com- 
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Fig.  24  (a)  Illustration  of  a  24  cell  5,h-order  quadrilateral  discretization  of  a  sphere,  (b)  Relative  mean  error  in  the 
RCS  of  a  sphere  of  radius  koa  =  6  computed  via  the  LCN  method  with  ninth-order  basis  versus  the 
interpolation  order  n  of  the  24  curvilinear  cells.  Also  compared  in  the  graph  is  the  relative  error  of  the  area 
of  the  patched  sphere. 


plex  geometries  with  a  high-order  method,  it  is  imperative  to  have  a  high-order  mesh  description. 

Currently  there  are  no  commercially  available  CAD  tools  that  provide  true  high-order 
modeling.  To  efficiently  model  complex  geometries  with  Mat-Scat,  a  high-order  mesh 
description  is  necessary.  Thus,  we  have  developed  a  mesh  processing  tool  that  can  generate  a 
high-order  mesh,  to  arbitrary  order,  from  a  low-order  mesh  generated  by  a  commercial  CAD 
tool.  The  mesh  processing  tool  is  named  UKY-MeshTool.  The  software  has  a  user  friendly  GUI 
interface  written  in  JAVA  (see  Fig.  25). 

The  UKY-MeshTool,  provides  high-order  surface/volume  meshes  for: 
o  Surfaces:  high-order  quadrilateral  or  triangular  elements, 
o  Volumes:  high-order  hexahedral,  prism,  or  tetrahedral  elements. 

It  is  emphasized  that  MeshTool  does  not  generate  the  initial  discretization.  Rather,  it  reads  in 
a  linear  mesh  (quadrilateral,  triangle,  hexahedral,  or  tetrahedral)  generated  by  another  CAD  tool 
Currently,  the  import  format  used  is  the  universal  file  format  of  SDRC  iDEAS™.  A  high-order 
mesh  is  then  generated  in  one  of  two  ways.  For  canonical  geometries,  where  individual  surfaces 
can  be  expressed  in  an  analytical  form,  a  scripting  software  can  be  used  to  provide  a  surface 
description  that  MeshTool  uses  to  generate  a  high-order  mesh  (c.fi,  Fig.  26).  The  second  manner 
is  to  read  in  a  coarse  and  fine-mesh  of  the  same  geometry  from  the  CAD  software.  This  requires 
the  CAD  software  to  support  a  “mesh-refinement”  feature  that  refines  individual  elements  of  a 
given  coarse  mesh.  MeshTool  does  not  require  any  specification  of  hierarchy  of  elements. 
Rather,  it  locates  nodes  shared  by  the  coarse  and  fine  meshes,  and  then  automatically  constructs 
the  hierarchy  of  elements  and  builds  the  higher-order  mesh. 
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Fig.  25 .  UKY-Mesh Tool  GUI  Interface  for  generating  and  processing  high-order  curvilinear  meshes. 


<2  ones. 
80  elements ' 


File  Eifet  Dtsjriay  ’  Functions  Checks 

Higher  Order  Elements  Build  From  Script, 

Higher  Order  {demerits  Build  From  O*  Script 
Higher  Order  Elements  Build  From  CoarseFine  Discretizations. 
Extrude  Surface 
Offset  Surface 

Extract  Surfaces  from  Selected  Volume  Cells 
Create  Finite  Periodic  Model 
FUp  Selected  Surfaces 

Enforce  selected  elernerrt  handedness  to  neighbors 
Edit  Material  Properties 
Change  Background  Properties  of  Elements 


Jsi*J 


7'">:  n  ■■■■ :  - 

m 


Fig.  26. Functions  supported  by  UKY-McshTool 
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;ig.  27. Mesh  checking  supported  by  UKY-McshTool. 
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Mesh  tool  also  allows  the  user  to  extrude  volume  elements  off  of  a  surface.  This  is  useful  for 
creating  coatings  on  a  surface,  extending  a  substrate  or  superstate,  etc.  Again,  curvilinear  offset 
elements  are  supported  via  a  scripted  description.  Alternatively,  surface  elements  can  be  offset 
in  a  similar  way  for  a  SIE  formulation.  Mesh  tool  will  also  allow  the  user  to  superimpose 
multiple  model  geometries.  This  allows  the  user  to  re-use  models,  or  to  facilitate  the  building  of 
more  complex  structures. 

Once  a  mesh  is  generated,  some  mesh  quality  checking  is  supported  by  MeshTool  such  as 
duplicate  node  checking  (c.f.,  Fig.  27).  It  will  also  compare  physical  parameters  with 
background  parameters  and  eliminate  any  volume  elements  where  these  are  identical  since  they 
will  support  a  zero  current.  Also,  when  modeling  surfaces,  the  direction  of  the  normal  vector  as 
determined  by  the  reciprocal  unitary  vector  a3  of  the  surface  patch  is  important  to  identify 
materials  above  and  below.  MeshTool  will  check  for  consistency  of  surface  normal  definition, 
and  can  flip  surface  normals  locally  or  globally  so  that  normals  are  consistent  with  the  surface 
parameterization. 

MeshTool  us  es  s  tandard  Legrange  int  erpolation  f  or  a  geometry  de  scription.  F  or  e  xample, 
aquadrilateral  patch  of  order  n  is  represented  by  («  +  l)x(«  +  l)  nodes  that  lie  on  the  surface.  In 

the  unitary  space,  these  nodes  are  uniformly  spaced,  as  illustrated  in  Fig.  28.  Each  node  has  a 
physical  coordinate  rjk  ( j,k  =  0,n ).  The  position  at  any  arbitrary  coordinate  ( u\u 2)  can  be 

obtained  via  interpolation: 

(*>:(•<%  <65) 

k-0  ./=0 

where  the  interpolation  polynomials  are  expressed  as: 

®"(u)  =  R'(n,u)R„_,(n,  1-w)  (66) 

where  Rfn,u)  is  a  Sylvester  interpolation  polynomial  [38]: 


R-,{n,u)  =  < 


1  •  k=  0 

1, 


1  <  /  <  n 
i  =  0 


(67) 


It  is  noted  that  this  interpolation  procedure  exactly  represents  a  bi-linear  quadrilateral  when 
n  =  1 ,  bi-quadratic  quadrilateral  when  n  =  2 ,  and  a  bi-cubic  quadrilateral  when  n  -  3 .  It  also 
represents  interpolations  to  arbitrary  order. 
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The  unitary  vectors  for  each  patch  are  also  readily  computed  as: 


du'  k=0  o 


(68) 


Analytical  expressions  for  the  derivatives  of  the  interpolation  polynomials  are  easily  derived,  and 
the  derivatives  of  the  Sylvester  interpolation  polynomials  can  be  expressed  via  a  recursive 
relationship  that  is  efficiently  computed. 

This  scheme  is  general  enough  to  interpolate  surfaces  up  to  arbitrary  order.  However,  it  still 
must  be  realized  that  the  interpolation  scheme  is  only  C°  continuous.  And,  for  very  high-orders 
(~ n  >  10)  the  interpolation  operator  can  be  ill  conditioned.  As  a  result,  Gibb’s  phenomena  can 
occur  leading  to  small  high  spatial-frequency  oscillation  of  the  interpolated  surface.  One  way 
that  this  has  been  addressed  is  that  MeshTool  can  generate  high-order  meshes  based  on  a  non- 
uniform  Gauss-Lobatto  point  spacing.  Nevertheless,  the  uniformly  spaced  interpolation  scheme 
proposed  in  (65)  works  quite  well  to  reasonably  high  order  (~  10th  order).  The  non-uniform 
spaced  meshing  has  shown  an  improvement  in  accuracy  for  higher-order  meshes. 


Finally,  MeshTool  exports  the  mesh  for  Mat-Scat  and  allows  the  user  to  define  source 
excitation  and  post-processing  directives  (c.f.,  Fig.  29). 
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5  Technology  Transfer 

5.A.  Publications 

Peer  Reviewed  Journals 

The  following  papers  published  in  peer  reviewed  journals  have  reported  in  the  open  literature  the 

results  and  methodologies  developed  under  this  contract: 

Dissertations 

[1]  S.  D.  Gedney,  "High-Order  Method  of  Moment  Solution  of  the  Scattering  by  Three- 
Dimensional  PEC  Bodies  using  Quadrature  Based  Point  Matching,"  Microwave  and 
Optical  Technology’  Letters,  vol.  29,  pp.  303-309,  June  5,  2001.  (5) 

[2]  Stephen  D.  Gedney,  Aiming  Zhu  ,  Wee-Hua  Tang,  Gang  Liu,  and  Peter  Petre  ,  “A  Fast, 
High-Order  Quadrature  Sampled  Pre-Corrected  FFT  for  Electromagnetic  Scattering,” 
Microwave  and  Optical  Technology  Letters,  vol.  36,  no.  5,  pp.  343-349,  March  5,  2003.  (2) 

[3]  S.  D.  Gedney  and  C.  C.  Lu,  “High-Order  Solution  for  the  Electromagnetic  Scattering  by 
Inhomogeneous  Dielectric  Bodies,”  Radio  Science,  vol.  38,  no.  1,  art.  no.  1015,  2003. 

[4]  G.  Liu  and  S.  D.  Gedney,  “High-Order  Moment  Method  Solution  for  the  Scattering 
Analysis  of  Penetrable  Bodies,”  Electromagnetics,  vol.  23,  no.  4,  pp.  331-346,  2003. 

[5]  S.  D.  Gedney,  "On  Deriving  a  Locally  Corrected  Nystrom  Scheme  from  a  Quadrature 
Sampled  Moment  Method,"  IEEE  Transactions  on  Antennas  and  Propagation,  vol.  51,  no. 
9,  pp.  2402-2412,  Sept.  2003. 

[6]  A.  Zhu  and  S.  D.  Gedney,  “A  Quadrature  Sampled  Pre-Corrected  FFT  for  the 
Electromagnetic  Scattering  from  Inhomogeneous  Objects,”  IEEE  Antennas  and  Wireless 
Propagation  Letters,  Vol.  2,  no.  1,  pp.  50-53,  2003. 

[7]  S.  D.  Gedney,  A.  Zhu,  and  C.  C.  Lu,  “Study  of  Mixed-Order  Basis  Functions  for  the 
Locally-Corrected  Nystrom  Method,”  IEEE  Transactions  on  Antennas  and  Propagation, 
vol.  52  no.  1,  pp.  (in  press),  November  2004. 

[8]  S.  D.  Gedney,  “Implementing  the  Locally  Corrected  Nystrom  method,”  Applied 
Computational  Electromagnetics  Society  Newsletter,  Vol.  18,  no.  3,  pp.  15-27,  November 
2003. 

Conference  Proceedings  -  Full  Paper 

The  following  papers  published  as  full  papers  in  peer  reviewed  conference  proceedings  have 

reported  in  the  open  literature  the  results  and  methodologies  developed  under  this  contract: 

[1]  S.  Gedney,  A.  Zhu,  W.  H.  Tang,  and  P.  Petre,  “High-Order  Pre-Corrected  FFT  Solution  for 
Electromagnetic  Scattering,”  2002  IEEE  International  Symposium  on  Antennas  and 
Propagation,  San  Antonio,  TX,  June  16-21,  2002. 

[2]  S.  Gedney  and  C.  C.  Lu,  “High-Order  Integral  Equation  Solution  for  Scattering  by 
Penetrable  Inhomogeneous  Volumes,”  2002  IEEE  International  Symposium  on  Antennas 
and  Propagation,  San  Antonio,  TX,  June  16-21,  2002. 

[3]  S.  Gedney  &  C.  C.  Lu,  “High-Order  Integral  Equation  Solution  Based  On  a  Hybrid 
Volume/Surface  Formulation,”  The  Annual  Review  of  Progress  in  Applied  Computational 
Electromagnetics,  Monterey,  CA,  March  24-28,  2003. 
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[4]  A.  Zhu  and  S.  Gedney,  ”  A  Fast,  High-Order  Integral  Equation  Solution  for  the  Scattering 
by  Inhomogeneous  Objects,”  2003  IEEE  International  Symposium  on  Antennas  and 
Propagation,  Columbus,  OH,  June  23-27 ,  vol.  1,  pp.  7-10,  2003. 

[5]  S.  D.  Gedney,  R.  Hannemann,  J.  Hannemann,  G.  Liu,  and  P.  Petre,  “A  Fast  Integral 
Equation  Solution  Technique  for  Printed  Circuits  in  Layered  Media,”  2003  IEEE 
International  Symposium  on  Antennas  and  Propagation,  Columbus,  OH,  June  23-27,  pp.  3 
-6  vol.  1,2003. 

[6]  S.  D.  Gedney  and  C.  C.  Lu,  “High-Order  Integral  Equation  Solution  for  Scattering  by 
Composite  Materials,”  2003  IEEE  International  Symposium  on  Antennas  and  Propagation, 
Columbus,  OH,  June  23-27,  pp.  1055  -  1058  vol.2,  2003. 

[7]  S.  D.  Gedney,  A.  Zhu,  C.  C.  Lu,  “High-order  Locally  Corrected  Nystrom  Solution  with 
Mixed-order  Basis  Functions  for  Electromagnetic  Scattering”,  ACES  Symposium,  Syracuse, 
NY,  April  19-23,  paper  no.  904178,  pg.  1-6,  2004 

[8]  A.  Zhu,  S.  D.  Gedney,  “Comparison  of  Muller  and  PMCHWT  Surface  Integral 
Formulations  for  the  Locally  Corrected  Nystrom  Method”,  2004  IEEE  International 
Symposium  on  Antennas  and  Propagation,  Monterey,  CA  June  21-25,  pp.  3871-3874,  2004 

[9]  A.  Zhu,  S.  D.  Gedney,  C.  Lu,  “Fast,  High-order,  Hybrid  Integral  Equation  Solver  for 
Electromagnetic  Scattering,”  2004  IEEE  International  Symposium  on  Antennas  and 
Propagation,  Monterey,  CA  June  21-25,  pp.  1 199  -  1202,  2004 

[10]  L.  Xuan,  A.  Zhu,  R.  J.  Adams,  S.  D.  Gedney,  “A  Broad  Band  Multilevel  Fast  Multipole 
Algorithm”,  2004  IEEE  International  Symposium  on  Antennas  and  Propagation, 
Monterey,  CA  June  21-25,  pp.  1195-1 198,  2004 

[11]  S.  D.  Gedney,  A.  Zhu,  C.  C.  Lu,  "Mixed-Order  Basis  Functions  for  the  Locally-Corrected 
Nystrom  M  ethod, "  2004  IEEE  International  Symposium  on  Antennas  and  Propagation, 
Monterey,  CA  June  21-25,  pp.  4044  -  4047,  2004 


The  following  dissertation  has  resulted  from  the  funding  of  this  program: 

Aiming  Zhu,  Ph.D.  Dissertation,  Fast,  Hybrid  High-Order  Integral  Solver  for  Electromagnetic 
Scattering,  June  2004. 


5.B.  Inventions 

No  patents  have  been  applied  for  from  this  research. 
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